import gdal
from gdal import ogr
from gdal import osr
import zipfile
import shutil
import json
from tqdm import tqdm
import re
import os
from basicTool import UTL


class ShapeFile:

    def __init__(self):
        pass

    @staticmethod
    def LayerIndexToLayerName(shapefilepath, layerindex=0):
        """
        根据图层索引获取图层名称
        :param shapefilepath:
        :param layerindex:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        if layerindex:
            layer = datasource.GetLayer(layerindex)
            layername = layer.GetName()
        else:
            layer = datasource.GetLayer(0)
            layername = layer.GetName()
        datasource.Destroy()
        return layername

    @staticmethod
    def GdalOgrDataTypeToPyhtonDataType(gdalogrdatatypename, gdalogrvalue):
        """
        gdal ogr 转为python的数据类型
        :param gdalogrdatatypename:
        :param gdalogrvalue:
        :return:
        """
        if gdalogrdatatypename == "Real":
            pythonvalue = float(gdalogrvalue)
        else:
            pythonvalue = int(gdalogrvalue)
        return pythonvalue

    @staticmethod
    def GdalOgrDataTypeNameToPythonDataTypeName(gdalogrdatatypename):
        """
        由gdal ogr的数据类型转换为python的数据类型名
        :param gdalogrdatatypename:
        :return:
        """
        if gdalogrdatatypename == "Real":
            pythondatatypename = "float"
        else:
            pythondatatypename = "int"
        return pythondatatypename

    @staticmethod
    def GetFieldDataType(shapefilepath, layername, fieldname):
        """
        获取shape file 的指定字段数据类型
        :param shapefilepath: 
        :param layername: 
        :param fieldname: 
        :return: 
        """
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        featuredefn = layer.GetLayerDefn()
        fieldindex = featuredefn.GetFieldIndex(fieldname)
        fielddefn = featuredefn.GetFieldDefn(fieldindex)
        datatype = fielddefn.GetType()
        datatypename = fielddefn.GetTypeName()
        datasource.Destroy()
        return datatypename

    @staticmethod
    def OpenGDBFile(gdbfilepath):
        # driver = ogr.GetDriverByName("File GDB")
        # gdbdriver = ogr.GetDriverByName("FileGDB")
        datasource = ogr.Open(gdbfilepath)
        driver = datasource.GetDriver()
        drivername = driver.GetName()
        print(drivername)
        layercount = datasource.GetLayerCount()
        for layerindex in range(layercount):
            layer = datasource.GetLayer(layerindex)
            layername = layer.GetName()
            print(layername)
        return

    @staticmethod
    def GetShapeFileSpatialreferencesystemWkt(shapefilepath, layername=None):
        """
        导出shape file 的空间参考系为proj 的json格式
        :param shapefilepath:
        :param layername:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        spatialreferencesystem = layer.GetSpatialRef()
        spatialreferencesystemwkt = spatialreferencesystem.ExportToWkt()
        datasource.Destroy()
        return spatialreferencesystemwkt

    @staticmethod
    def GetShapeFileSpatialreferencesystemProjJson(shapefilepath, layername=None):
        """
        导出shape file 的空间参考系为proj 的json格式
        :param shapefilepath:
        :param layername:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        spatialreferencesystem = layer.GetSpatialRef()
        spatialreferencesystemprojjson = spatialreferencesystem.ExportToPROJJSON()
        datasource.Destroy()
        return spatialreferencesystemprojjson

    @staticmethod
    def GetShapeFileSpatialreferencesystemProj4(shapefilepath, layername=None):
        """
        导出shape file 的空间参考系为proj 的json格式
        :param shapefilepath:
        :param layername:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        spatialreferencesystem = layer.GetSpatialRef()
        spatialreferencesystemproj4 = spatialreferencesystem.ExportToProj4()
        datasource.Destroy()
        return spatialreferencesystemproj4
    
    @staticmethod
    def GetShapeFileLayerCount(shapefilepath):
        """
        获取shape file 的layer 数量
        :param shapefilepath:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        layercount = datasource.GetLayerCount()
        datasource.Destroy()
        return layercount

    @staticmethod
    def GetShapefileFeatureDictList(shapefilepath, layername=None):
        """
        获得shp file 要素字典列表
        :param shapefilepath:
        :param layername:
        :return:
        """
        gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        featurecount = layer.GetFeatureCount()
        featuredictlist = []
        for featureindex in range(featurecount):
            feature = layer.GetFeature(featureindex)
            featurejson = feature.ExportToJson()
            featuredict = json.loads(featurejson)
            featuredictlist.append(featuredict)
            feature.Destroy()
        datasource.Destroy()
        return featuredictlist

    @staticmethod
    def GetShapeFileLayerNameList(shapefilepath):
        """
        获取shape file 的 layername list表
        :param shapefilepath:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        layercount = datasource.GetLayerCount()
        layernamelist = []
        for layerindex in range(layercount):
            layer = datasource.GetLayer(layerindex)
            layername = layer.GetName()
            layernamelist.append(layername)
        datasource.Destroy()
        return layernamelist

    @staticmethod
    def CreateEmptyShapeFile(shapefilepath, layername, spatialreferencesystemprojjson, geometrytype):
        """
        生成空的shp file
        :param shapefilepath:
        :param layername:
        :param spatialreferencesystemprojjson:
        :param geometrytype:
        :return:
        """
        drivername = "ESRI Shapefile"
        driver = ogr.GetDriverByName(drivername)
        datasource = driver.CreateDataSource(shapefilepath)
        spatialreferencesystem = osr.SpatialReference()
        spatialreferencesystem.ImportFromProj4(spatialreferencesystemprojjson)
        layer = datasource.CreateLayer(layername, srs=spatialreferencesystem, geom_type=geometrytype)
        datasource.Destroy()
        return
    
    @staticmethod
    def GetShapeFileGeometrytype(shapefilepath, layername=None):
        """
        获得shp file 的几何类型
        :param shapefilepath:
        :param layername:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        geometrytype = layer.GetGeomType()
        datasource.Destroy()
        return geometrytype

    @staticmethod
    def GetShapefileFieldList(shapefilepath, layername=None):
        """
        获取shp file 字段定义列表
        :param shapefilepath:
        :param layername:
        :return:
        """
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        featuredefn = layer.GetLayerDefn()
        fieldcount = featuredefn.GetFieldCount()
        fielddictlist = []
        for fieldindex in range(fieldcount):
            fielddefn = featuredefn.GetFieldDefn(fieldindex)
            fielddict = {
                "name": fielddefn.name,
                "type": fielddefn.type,
                "width": fielddefn.width,
                "justify": fielddefn.justify,
                "precision": fielddefn.precision}
            fielddictlist.append(fielddict)
        datasource.Destroy()
        return fielddictlist

    @staticmethod
    def GetSpatialReferenceSystemType(shapefilepath, layername=None):
        datasource = ogr.Open(shapefilepath)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        spatialreferencesystem = layer.GetSpatialRef()
        if spatialreferencesystem.IsGeographic():
            spatialreferencesystemname = spatialreferencesystem.GetName()
            print("投影坐标系", "坐标系名称为：", spatialreferencesystemname)
        elif spatialreferencesystem.IsProjected():
            spatialreferencesystemname = spatialreferencesystem.GetName()
            print("地理坐标系", "坐标系名称为：", spatialreferencesystemname)
        else:
            print("非坐标系")
        datasource.Destroy()
        return

    @staticmethod
    def FeatureDictListTransformSpatialreferenceSystem(
            featuredictlist, originalspatialreferencesystemwkt, targetspatialreferencesystemwkt):
        """
        转换要素列表几何坐标（投影转地理）
        :param featuredictlist:
        :param originalspatialreferencesystemwkt:
        :param targetspatialreferencesystemwkt:
        :return:
        """
        originalspatialreferencesystem = osr.SpatialReference()
        originalspatialreferencesystem.ImportFromWkt(originalspatialreferencesystemwkt)
        targetspatialreferencesystem = osr.SpatialReference()
        targetspatialreferencesystem.ImportFromWkt(targetspatialreferencesystemwkt)
        coordinatetransformation = osr.CoordinateTransformation(
            originalspatialreferencesystem, targetspatialreferencesystem)
        newfeaturedictlist = []
        for featuredict in featuredictlist:
            newfeaturedict = featuredict
            geometrydict = featuredict["geometry"]
            geometrytype = geometrydict["type"]
            if geometrytype == "Polygon":
                polygon = geometrydict["coordinates"]
                newpolygon = []
                for ring in polygon:
                    newring = []
                    for point in ring:
                        pointtuple = coordinatetransformation.TransformPoint(point[0], point[1])
                        newpoint = [pointtuple[1], pointtuple[0]]
                        newring.append(newpoint)
                    newpolygon.append(newring)
                newfeaturedict["geometry"]["coordinates"] = newpolygon
                newfeaturedictlist.append(newfeaturedict)
            elif geometrytype == "MultiPolygon":
                fielddictlist = featuredict["properties"]
                polygonlist = geometrydict["coordinates"]
                newpolygonlist = []
                for polygon in polygonlist:
                    newpolygon = []
                    for ring in polygon:
                        newring = []
                        for point in ring:
                            pointtuple = coordinatetransformation.TransformPoint(point[0], point[1])
                            newpoint = [pointtuple[1], pointtuple[0]]
                            newring.append(newpoint)
                        newpolygon.append(newring)
                    newpolygonlist.append(newpolygon)
                newfeaturedict["geometry"]["coordinates"] = newpolygonlist
                newfeaturedictlist.append(newfeaturedict)
            elif geometrytype == "Point":
                print(geometrytype)
            else:
                print("error!")
        return newfeaturedictlist

    @staticmethod
    def FeatureDictListTransformSpatialreferenceSystem2(
            featuredictlist, originalspatialreferencesystemwkt, targetspatialreferencesystemwkt):
        """
        转换要素列表几何坐标（地理转投影）
        :param featuredictlist:
        :param originalspatialreferencesystemwkt:
        :param targetspatialreferencesystemwkt:
        :return:
        """
        originalspatialreferencesystem = osr.SpatialReference()
        originalspatialreferencesystem.ImportFromWkt(originalspatialreferencesystemwkt)
        targetspatialreferencesystem = osr.SpatialReference()
        targetspatialreferencesystem.ImportFromWkt(targetspatialreferencesystemwkt)
        coordinatetransformation = osr.CoordinateTransformation(
            originalspatialreferencesystem, targetspatialreferencesystem)
        newfeaturedictlist = []
        for featuredict in featuredictlist:
            newfeaturedict = featuredict
            geometrydict = featuredict["geometry"]
            geometrytype = geometrydict["type"]
            if geometrytype == "Polygon":
                polygon = geometrydict["coordinates"]
                newpolygon = []
                for ring in polygon:
                    newring = []
                    for point in ring:
                        pointtuple = coordinatetransformation.TransformPoint(point[1], point[0])
                        newpoint = [pointtuple[0], pointtuple[1]]
                        newring.append(newpoint)
                    newpolygon.append(newring)
                newfeaturedict["geometry"]["coordinates"] = newpolygon
                newfeaturedictlist.append(newfeaturedict)
            elif geometrytype == "MultiPolygon":
                fielddictlist = featuredict["properties"]
                polygonlist = geometrydict["coordinates"]
                newpolygonlist = []
                for polygon in polygonlist:
                    newpolygon = []
                    for ring in polygon:
                        newring = []
                        for point in ring:
                            pointtuple = coordinatetransformation.TransformPoint(point[1], point[0])
                            newpoint = [pointtuple[0], pointtuple[1]]
                            newring.append(newpoint)
                        newpolygon.append(newring)
                    newpolygonlist.append(newpolygon)
                newfeaturedict["geometry"]["coordinates"] = newpolygonlist
                newfeaturedictlist.append(newfeaturedict)
            elif geometrytype == "Point":
                print(geometrytype)
            else:
                print("error!")
        return newfeaturedictlist

    @staticmethod
    def FieldDictListUpdateShapeFile(fielddictlist, shapefilepath, layername=None):
        """
        依据
        :param fielddictlist:
        :param shapefilepath: 
        :param layername: 
        :return: 
        """
        datasource = ogr.Open(shapefilepath, 1)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        for fielddict in fielddictlist:
            field = ogr.FieldDefn()
            field.name = fielddict["name"]
            field.type = fielddict["type"]
            field.width = fielddict["width"]
            field.justify = fielddict["justify"]
            field.precision = fielddict["precision"]
            layer.CreateField(field)
            field.Destroy()
        datasource.Destroy()
        return
    
    @staticmethod
    def FeatureDictListUpdateShapeFile(featuredictlist, shapefilepath, layername=None):
        """
        根据要素字典列表更新shp file
        :param featuredictlist: 
        :param shapefilepath: 
        :param layername: 
        :return: 
        """
        # gdal.SetConfigOption("SHAPE_ENCODING", "UTF-8")
        gdal.SetConfigOption("SHAPE_ENCODING", "GBK")
        # gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES")
        datasource = ogr.Open(shapefilepath, 1)
        if layername:
            layer = datasource.GetLayerByName(layername)
        else:
            layer = datasource.GetLayer(0)
        layerdefn = layer.GetLayerDefn()
        fieldcount = layerdefn.GetFieldCount() 
        for featuredict in featuredictlist:
            feature = ogr.Feature(layerdefn)
            fielddict = featuredict["properties"]
            for key, value in fielddict.items():
                feature.SetField(key, value)
            geometrydict = featuredict["geometry"]
            geometryjson = json.dumps(geometrydict)
            geometry = ogr.CreateGeometryFromJson(geometryjson)
            feature.SetGeometry(geometry)
            layer.CreateFeature(feature)
        datasource.Destroy()
        return

    @staticmethod
    def LayerIntersect(firstshapefilepath, firstlayername, secondshapefilepath, secondlayername):
        """
        求两个图层中相交的要素
        :param firstshapefilepath:
        :param firstlayername:
        :param secondshapefilepath:
        :param secondlayername:
        :return:
        """
        firstdatasource = ogr.Open(firstshapefilepath)
        firstlayer = firstdatasource.GetLayerByName(firstlayername)
        firstfeaturecount = firstlayer.GetFeatureCount()
        seconddatasource = ogr.Open(secondshapefilepath)
        secondlayer = seconddatasource.GetLayerByName(secondlayername)
        secondfeaturecount = secondlayer.GetFeatureCount()
        countysentinel2featuredictlist = []
        for firstfeatureindex in range(firstfeaturecount):
            firstfeature = firstlayer.GetFeature(firstfeatureindex)
            firstgeometry = firstfeature.geometry()
            for secondfeatureindex in range(secondfeaturecount):
                secondfeature = secondlayer.GetFeature(secondfeatureindex)
                secondgeometry = secondfeature.geometry()
                if firstgeometry.Intersect(secondgeometry):
                    firstfeaturejson = firstfeature.ExportToJson()
                    firstfeaturedict = json.loads(firstfeaturejson)
                    countysentinel2featuredictlist.append(firstfeaturedict)
                    break
                else:
                    pass
        firstdatasource.Destroy()
        seconddatasource.Destroy()
        return countysentinel2featuredictlist


class RasterFile:

    def __init__(self):
        pass

    @staticmethod
    def GetRasterFileDataType(rasterfilepath):
        """
        获取栅格波段数据类型，并返GDAL数据类型对象
        :param rasterfilepath:
        :return:
        """
        raster = gdal.Open(rasterfilepath, 0)
        firstband = raster.GetRasterBand(1)
        datatype = firstband.DataType
        del raster
        return datatype

    @staticmethod
    def GetRasterFileNoDataValue(rasterfilepath):
        raster = gdal.Open(rasterfilepath, 0)
        firstband = raster.GetRasterBand(1)
        nodatavalue = firstband.GetNoDataValue()
        del raster
        return nodatavalue

    @staticmethod
    def GetRasterFileSpatialreferencesystemwkt(rasterfilepath):
        raster = gdal.Open(rasterfilepath)
        spatialreferencesystem = raster.GetSpatialRef()
        spatialreferencesystemwkt = spatialreferencesystem.ExportToWkt()
        del raster
        return spatialreferencesystemwkt

    @staticmethod
    def GetRasterFileGeoTransform(rasterfilepath):
        """
        获取栅格仿射参数，并返回元组
        :param rasterfilepath:
        :return:
        """
        raster = gdal.Open(rasterfilepath, 0)
        geotransform = raster.GetGeoTransform()
        del raster
        return geotransform

    @staticmethod
    def GetRasterFileProjection(rasterfilepath):
        """
        获取栅格投影信息，并返回字符串
        :param rasterfilepath:
        :return:
        """
        raster = gdal.Open(rasterfilepath, 0)
        project = raster.GetProjection()
        del raster
        return project

    @staticmethod
    def GetRasterFileShape(rasterfilepath):
        """
        获取栅格大小
        :param rasterfilepath:
        :return:
        """
        raster = gdal.Open(rasterfilepath)
        width, height, bandcount = raster.RasterXSize, raster.RasterYSize, raster.RasterCount
        shape = [width, height, bandcount]
        del raster
        return shape

    @staticmethod
    def PythonDataTypeNameToGdalDataTypeName(pythondatatypename):
        """
        由python的数据类型的名称获得gdal数据类型名称
        :param pythondatatypename: 
        :return: 
        """
        if pythondatatypename == "int":
            gdaldatatypename = "UInt16"
        elif pythondatatypename == "float":
            gdaldatatypename = "Float32"
        else:
            gdaldatatypename = "Byte"
        return gdaldatatypename
    
    @staticmethod
    def PythonDataTypeNameToGdalDataType(pythondatatypename):
        """
        由python数据类型名称获得gdal数据类型对象
        :param pythondatatypename: 
        :return: 
        """
        if pythondatatypename == "int":
            return gdal.GDT_UInt16
        elif pythondatatypename == "float":
            return gdal.GDT_Float32
        else:
            return gdal.GDT_Unknown
        
    @staticmethod
    def CreateSingleBandRasterFile(
            rasterfilepath, bandarray, spatialreferencesystemwkt, datetype):
        driver = gdal.GetDriverByName("GTiff")
        (channel, height, width) = bandarray.shape()
        spatialreferencesystem = osr.SpatialReference()
        spatialreferencesystem.ImportFromWkt(spatialreferencesystemwkt)
        dataset = driver.Create(rasterfilepath, width, height, 2, datetype)
        dataset.SetGeoTransform(geotrs)
        dataset.SetProjection(proj)
        dataset.SetSpatialRef(spatialreferencesystem)

        bandpositive = dataset.GetRasterBand(1)
        bandpositive.SetNoDataValue(0)
        bandpositive.FlushCache()
        gdal.RasterizeLayer(dataset, [1], layer, options=["ATTRIBUTE=BL"])
        gdal.RasterizeLayer(dataset, [1], layerpositive, burn_values=[1])
        gdal.RasterizeLayer(dataset, [2], layernegative, burn_values=[1])
        del dataset

        datasource.DeleteLayer(self.layerpositivename)
        datasource.DeleteLayer(self.layernegativename)
        datasource.Destroy()
        return


class Sentinel2:

    def __init__(self):
        pass

    @classmethod
    def Sentinel2ZipfileToRasterfile(cls, sentinel2zipfilepath, unzipdircetory):
        """
        由哨兵2压缩文件直接生成多波段栅格文件
        :param sentinel2zipfilepath: 哨兵2压缩文件
        :param unzipdircetory: 解压文件目录
        :return:
        """
        multiplesinglesentinel2rasterfilefolder = cls.Sentinel2ZipfileToUnzipfolder(sentinel2zipfilepath,
                                                                                    unzipdircetory)
        cls.Zone8BandResampling(multiplesinglesentinel2rasterfilefolder)
        zonerasterfilepath = cls.Zone8BandCompose(multiplesinglesentinel2rasterfilefolder, unzipdircetory)
        shutil.rmtree(multiplesinglesentinel2rasterfilefolder)
        return zonerasterfilepath
    
    @classmethod
    def Sentinel2SingleBandRaster(cls, sentinel2zipfilepath, unzipdircetory):
        """
        获取单波段影像
        :param sentinel2zipfilepath: 
        :param unzipdircetory: 
        :return: 
        """
        multiplesinglesentinel2rasterfilefolder = cls.Sentinel2ZipfileToUnzipfolder(
            sentinel2zipfilepath, unzipdircetory)
        cls.Zone8BandResampling(multiplesinglesentinel2rasterfilefolder)
        rasterfilenamelist = os.listdir(multiplesinglesentinel2rasterfilefolder)
        for rasterfilename in rasterfilenamelist:
            print(rasterfilename)
            [filename, extend] = os.path.splitext(rasterfilename)
            print(filename)
            filenamestrlist = filename.split("_")
            newfilenamestrlist = filenamestrlist[2:] + [filenamestrlist[1][3:]]
            newfilename = "_".join(newfilenamestrlist)
            newrasterfilename = newfilename + extend
            rasterfilepath = os.path.join(multiplesinglesentinel2rasterfilefolder, rasterfilename)
            newrasterfilepath = os.path.join(multiplesinglesentinel2rasterfilefolder, newrasterfilename)
            os.rename(rasterfilepath, newrasterfilepath)
        return 
    
    @classmethod
    def BatchSentinel2RemoveRasterToSaveDirectory(cls, sentinel2zipdirectory, savedirectory):
        """
        批量解压哨兵2压缩文件为但波段栅格文件到指定目录
        :param sentinel2zipdirectory: 
        :param savedirectory: 
        :return: 
        """
        allsubdirectorynamelist = os.listdir(sentinel2zipdirectory)
        pattern = re.compile(r".zip")
        sentinel2zipfilepathlist = []
        for subdirectoryname in allsubdirectorynamelist:
            if pattern.search(subdirectoryname):
                sentinel2zipfilepath = os.path.join(sentinel2zipdirectory, subdirectoryname)
                sentinel2zipfilepathlist.append(sentinel2zipfilepath)
        for sentinel2zipfilepath in tqdm(sentinel2zipfilepathlist):
            cls.Sentinel2RemoveRasterToSaveDirectory(sentinel2zipfilepath, savedirectory)
        return
    
    @classmethod
    def Sentinel2RemoveRasterToSaveDirectory(cls, sentinel2zipfilepath, unzipdircetory):
        """
        获取单波段影像到指定文件目录
        :param sentinel2zipfilepath: 
        :param unzipdircetory: 
        :return: 
        """
        multiplesinglesentinel2rasterfilefolder = cls.Sentinel2ZipfileToUnzipfolder(
            sentinel2zipfilepath, unzipdircetory)
        cls.Zone8BandResampling(multiplesinglesentinel2rasterfilefolder)
        rasterfilenamelist = os.listdir(multiplesinglesentinel2rasterfilefolder)
        for rasterfilename in rasterfilenamelist:
            [filename, extend] = os.path.splitext(rasterfilename)
            filenamestrlist = filename.split("_")
            newfilenamestrlist = filenamestrlist[2:] + [filenamestrlist[1][3:]]
            newfilename = "_".join(newfilenamestrlist)
            newrasterfilename = newfilename + extend
            rasterfilepath = os.path.join(multiplesinglesentinel2rasterfilefolder, rasterfilename)
            newrasterfilepath = os.path.join(unzipdircetory, newrasterfilename)
            shutil.move(rasterfilepath, newrasterfilepath)
        shutil.rmtree(multiplesinglesentinel2rasterfilefolder)
        return

    @classmethod
    def Sentinel2ZipfileToUnzipfolder(cls, sentinel2zipfilepath, unzipdir):
        """解压哨兵2栅格文件按固定格式到至指定文件目录中"""
        '''获取哨兵2原始压缩文件名中的信息'''
        sentinel2rasterfileinfolist = cls.Sentinel2GetZipFileDate(sentinel2zipfilepath)
        '''构建解压文件夹'''
        unzipfolder = cls.Sentinel2BuildUnzipFolderPath(sentinel2zipfilepath, sentinel2rasterfileinfolist, unzipdir)
        '''获取压缩文件名列表'''
        zipfilenamelist = cls.Sentinel2UnzipRasterNames(sentinel2zipfilepath)
        # [print(singlezipfilename) for singlezipfilename in zipfilenamelist]
        '''解压列表中的文件'''
        cls.UnzipFileToFolder(sentinel2zipfilepath, unzipfolder, zipfilenamelist)
        '''文件重命名'''
        cls.Sentinel2RasterFilesRename(unzipfolder)
        return unzipfolder

    @staticmethod
    def Sentinel2GetZipFileDate(sentinel2zipfilepath):
        """
        获取哨兵2原始压缩文件名中的信息
        :param sentinel2zipfilepath: 哨兵2压缩文件
        :return:
        """
        if not os.path.exists(sentinel2zipfilepath):
            return print(sentinel2zipfilepath, "\n该路径不存在！")
        [zipfiledirectory, zipfilename] = os.path.split(sentinel2zipfilepath)
        [filename, extend] = os.path.splitext(zipfilename)
        filesignstrlist = filename.split("_")
        satellite, product = filesignstrlist[0], filesignstrlist[1]
        zoneid, date = filesignstrlist[-2], filesignstrlist[-1][:8]
        return [satellite, product, zoneid, date]

    @staticmethod
    def Sentinel2BuildUnzipFolderPath(sentinel2rasterfilezippath, sentinel2rasterfileinfolist, unzipdir):
        """
        依据指定格式创建哨兵2解压文件夹
        :param sentinel2rasterfilezippath:
        :param sentinel2rasterfileinfolist:
        :param unzipdir:
        :return:
        """
        [satellite, product, zoneid, date] = sentinel2rasterfileinfolist
        # year, month, day = int(date[:4]), int(date[4:6]), int(date[6:])
        # dayc = UTL.YearMonthDayToNumber(year, month, day)
        # foldername = "%s_%s_%s_%s_%d" % (satellite, product, zoneid, date, dayc)
        foldername = "%s_%s_%s_%s" % (satellite, product, zoneid, date)
        unzipfolder = os.path.join(unzipdir, foldername)
        if os.path.exists(unzipfolder):
            return print(unzipfolder, ":解压文件夹存在！")
        else:
            os.makedirs(unzipfolder)
        return unzipfolder

    @staticmethod
    def Sentinel2UnzipRasterNames(zipfilepath):
        """
        获取需要解压的哨兵2压缩文件中文件列表
        :param zipfilepath:
        :return:
        """
        allzipfile = zipfile.ZipFile(zipfilepath, "r")
        allfiles = allzipfile.namelist()
        '''获取压缩文件中的指定波段文件列表'''
        band8jpgfiles = []
        for i in range(1, 9):
            pattern = re.compile("_B0[%d](.)*\.jp2" % i)
            bandjp2s = [file for file in allfiles if re.search(pattern, file)]
            if len(bandjp2s) >= 2:
                rasfilepres = []
                for bandjp2 in bandjp2s:
                    jp2filename = os.path.split(bandjp2)[1]
                    filename = os.path.splitext(jp2filename)[0]
                    rasfilepres.append(int(filename.split("_")[-1][:-1]))
                goalindex = rasfilepres.index(min(rasfilepres))
                band8jpgfiles.append(bandjp2s[goalindex])
            else:
                band8jpgfiles.append(bandjp2s[0])
        return band8jpgfiles

    @staticmethod
    def UnzipFileToFolder(zipfilepath, unzipfolder, zipfilenames=None):
        """
        依据列表解压压缩文件中的文件到指定的新建文件夹下
        :param zipfilepath:
        :param unzipfolder:
        :param zipfilenames:
        :return:
        """
        allzipfile = zipfile.ZipFile(zipfilepath, "r")
        if zipfilenames is None:
            zipfilenames = allzipfile.namelist()
        for zipfilename in zipfilenames:
            allzipfile.extract(zipfilename, unzipfolder)
            file = os.path.join(unzipfolder, zipfilename)
            filename = os.path.split(zipfilename)[1]
            newfile = os.path.join(unzipfolder, filename)
            shutil.move(file, newfile)
        allzipfile.close()
        foldernames = os.listdir(unzipfolder)
        for foldername in foldernames:
            folder = os.path.join(unzipfolder, foldername)
            if os.path.isdir(folder):
                shutil.rmtree(folder)
        return

    @staticmethod
    def Sentinel2RasterFilesRename(unzipfolder):
        """
        依据文件夹的名字修改每个波段栅格文件名
        :param unzipfolder:
        :return:
        """
        filenames = os.listdir(unzipfolder)
        foldername = os.path.split(unzipfolder)[1]
        for i in range(1, 9):
            pattern = re.compile('_(.)*_B0%d(.)*\.jp2' % i)
            filenames = os.listdir(unzipfolder)
            for filename in filenames:
                if re.search(pattern, filename):
                    extend = os.path.splitext(filename)[1]
                    filepath = os.path.join(unzipfolder, filename)
                    newfilename = "%s_B0%d%s" % (foldername, i, extend)
                    newfilepath = os.path.join(unzipfolder, newfilename)
                    os.rename(filepath, newfilepath)
        return

    @staticmethod
    def Zone8BandResampling(zonesentinel2rasterfolder):
        """
        对提取后原始哨兵影像进行重采样，使所有栅格分辨率都与最高的栅格一致！
        :param zonesentinel2rasterfolder:
        :return:
        """
        sentinel2singlebandrasterfilepathlist = UTL.GetGoalFile(zonesentinel2rasterfolder, r".jp2")
        rasterfilenamelist = os.listdir(zonesentinel2rasterfolder)
        pres, geots = [], []
        for rasterfilename in rasterfilenamelist:
            rasterfilepath = os.path.join(zonesentinel2rasterfolder, rasterfilename)
            rasterfile = gdal.Open(rasterfilepath)
            geotransform = rasterfile.GetGeoTransform()
            pres.append(abs(geotransform[1]))
            geots.append(geotransform)
            del rasterfile
        standardindex = pres.index(min(pres))
        goalgeot = geots[standardindex]
        for sentinel2singlebandrasterfilepath in sentinel2singlebandrasterfilepathlist:
            UTL.RasterReshapeByGeo(sentinel2singlebandrasterfilepath, goalgeot, zonesentinel2rasterfolder)
            os.remove(sentinel2singlebandrasterfilepath)
        return

    @staticmethod
    def Zone8BandCompose(zonesentinel2rasterfolder, savefolder):
        """
        组合哨兵前8波段
        :param zonesentinel2rasterfolder:
        :param savefolder:
        :return:
        """
        singlerasterfilenamelist = os.listdir(zonesentinel2rasterfolder)
        bandcount = len(singlerasterfilenamelist)
        rasterfilenamesign = os.path.splitext(singlerasterfilenamelist[0])[0]
        signstrs = rasterfilenamesign.split("_")[:-1]
        zonerasterfilename = "_".join(signstrs)
        zonerastertiffilename = zonerasterfilename + ".tif"
        singlerasterfilepath = os.path.join(zonesentinel2rasterfolder, singlerasterfilenamelist[0])
        singlerasterfile = gdal.Open(singlerasterfilepath)
        width, height = singlerasterfile.RasterXSize, singlerasterfile.RasterYSize
        geotransform = singlerasterfile.GetGeoTransform()
        spatialcoordinatesystem = singlerasterfile.GetSpatialRef()
        projection = singlerasterfile.GetProjectionRef()
        datatype = singlerasterfile.GetRasterBand(1).DataType
        del singlerasterfile
        arraylist = []
        for index in range(1, bandcount + 1):
            singlerasterfilename = zonerasterfilename + "_B0%d.tif" % index
            singlerasterfilepath = os.path.join(zonesentinel2rasterfolder, singlerasterfilename)
            singlerasterfile = gdal.Open(singlerasterfilepath)
            arraylist.append(singlerasterfile.GetRasterBand(1).ReadAsArray(0, 0, width, height))
            del singlerasterfile
        zonerasterfilepath = os.path.join(savefolder, zonerastertiffilename)
        driver = gdal.GetDriverByName("GTiff")
        zonerasterfile = driver.Create(zonerasterfilepath, width, height, bandcount, datatype)
        zonerasterfile.SetGeoTransform(geotransform)
        zonerasterfile.SetSpatialRef(spatialcoordinatesystem)
        zonerasterfile.SetProjection(projection)
        [zonerasterfile.GetRasterBand(i + 1).WriteArray(arraylist[i]) for i in range(bandcount)]
        del zonerasterfile
        return zonerasterfilepath


def ShapeFileToRasterFile(newrasterfilepath, rasterfilepath, shapefilepath, layername=None):
    """
    矢量转栅格
    :param newrasterfilepath:
    :param rasterfilepath:
    :param shapefilepath:
    :param layername:
    :return:
    """
    driver = gdal.GetDriverByName("GTiff")
    bandcount = 1
    shape = RasterFile.GetRasterFileShape(rasterfilepath)
    datatype = RasterFile.GetRasterFileDataType(rasterfilepath)
    geotransform = RasterFile.GetRasterFileGeoTransform(rasterfilepath)
    project = RasterFile.GetRasterFileProjection(rasterfilepath)
    spatialreferencesystemwkt = RasterFile.GetRasterFileSpatialreferencesystemwkt(rasterfilepath)
    spatialreferencesystem = osr.SpatialReference()
    nodatavalue = RasterFile.GetRasterFileNoDataValue(rasterfilepath)
    spatialreferencesystem.ImportFromWkt(spatialreferencesystemwkt)
    [width, height] = shape[:2]
    dataset = driver.Create(newrasterfilepath, width, height, bandcount, datatype)
    dataset.SetGeoTransform(geotransform)
    dataset.SetProjection(project)
    dataset.SetSpatialRef(spatialreferencesystem)
    band = dataset.GetRasterBand(1)
    if nodatavalue is not None:
        band.SetNoDataValue(nodatavalue)
    band.FlushCache()
    datasource = ogr.Open(shapefilepath)
    if layername:
        layer = datasource.GetLayerByName(layername)
    else:
        layer = datasource.GetLayer(0)
    gdal.RasterizeLayer(dataset, [1], layer, burn_values=[1])
    datasource.Destroy()
    del dataset
    return


def ShapeToRasterFillFieldValue(
        newrasterfilepath, rasterfilepath, shapefilepath, layername, fieldname):
    """
    矢量转栅格
    :param newrasterfilepath:
    :param rasterfilepath:
    :param shapefilepath:
    :param layername:
    :param fieldname:
    :return:
    """
    driver = gdal.GetDriverByName("GTiff")
    bandcount = 1
    shape = RasterFile.GetRasterFileShape(rasterfilepath)
    datatype = RasterFile.GetRasterFileDataType(rasterfilepath)
    ogrdatatypename = ShapeFile.GetFieldDataType(shapefilepath, layername, fieldname)
    pythondatatypename = ShapeFile.GdalOgrDataTypeNameToPythonDataTypeName(ogrdatatypename)
    fieldgdaltype = RasterFile.PythonDataTypeNameToGdalDataType(pythondatatypename)
    geotransform = RasterFile.GetRasterFileGeoTransform(rasterfilepath)
    project = RasterFile.GetRasterFileProjection(rasterfilepath)
    spatialreferencesystemwkt = RasterFile.GetRasterFileSpatialreferencesystemwkt(rasterfilepath)
    spatialreferencesystem = osr.SpatialReference()
    nodatavalue = RasterFile.GetRasterFileNoDataValue(rasterfilepath)
    spatialreferencesystem.ImportFromWkt(spatialreferencesystemwkt)
    [width, height] = shape[:2]
    dataset = driver.Create(newrasterfilepath, width, height, bandcount, fieldgdaltype)
    dataset.SetGeoTransform(geotransform)
    dataset.SetProjection(project)
    dataset.SetSpatialRef(spatialreferencesystem)
    band = dataset.GetRasterBand(1)
    if nodatavalue is not None:
        band.SetNoDataValue(nodatavalue)
    band.FlushCache()
    datasource = ogr.Open(shapefilepath)
    if layername:
        layer = datasource.GetLayerByName(layername)
    else:
        layer = datasource.GetLayer(0)
    gdal.RasterizeLayer(dataset, [1], layer, options=["ATTRIBUTE=%s" % fieldname])
    datasource.Destroy()
    del dataset
    return


def ShapeFileTransformSpatialreferencesystem(
        newspatialreferencesystemwkt, shapefilepath, layername, newshapefilepath, newlayername):
    """
    shape file 空间参考系转换（地理转投影）
    :param newspatialreferencesystemwkt: 
    :param shapefilepath: 
    :param layername: 
    :param newshapefilepath: 
    :param newlayername: 
    :return: 
    """
    shapespatialreferencesystemwkt = ShapeFile.GetShapeFileSpatialreferencesystemWkt(
        shapefilepath, layername)
    geometrytype = ShapeFile.GetShapeFileGeometrytype(shapefilepath, layername)
    ShapeFile.CreateEmptyShapeFile(
        newshapefilepath, newlayername, newspatialreferencesystemwkt, geometrytype)
    fielddictlist = ShapeFile.GetShapefileFieldList(shapefilepath, layername)
    featuredictlist = ShapeFile.GetShapefileFeatureDictList(shapefilepath, layername)
    ShapeFile.FieldDictListUpdateShapeFile(fielddictlist, newshapefilepath, newlayername)
    newfeaturedictlist = ShapeFile.FeatureDictListTransformSpatialreferenceSystem2(
        featuredictlist, shapespatialreferencesystemwkt, newspatialreferencesystemwkt)
    ShapeFile.FeatureDictListUpdateShapeFile(newfeaturedictlist, newshapefilepath, newlayername)
    return


class ChinaSentinel2:

    def __init__(self):
        pass

    @staticmethod
    def TransFormWorldSentinel2ZoneShapeFile(
            worldsentinel2zoneshapefilepath, spatialreferencesystemshapefilepath, newshapefilepath):
        """
        根据全世界哨兵投影图幅shape file与wgs84的shape file生成地理单图层shape file 
        :param worldsentinel2zoneshapefilepath: 
        :param spatialreferencesystemshapefilepath: 
        :param newshapefilepath:
        :return: 
        """
        spatialreferencesystemwkt = ShapeFile.GetShapeFileSpatialreferencesystemWkt(spatialreferencesystemshapefilepath)
        layernamelist = ShapeFile.GetShapeFileLayerNameList(worldsentinel2zoneshapefilepath)
        ShapeFile.GetSpatialReferenceSystemType(spatialreferencesystemshapefilepath)
        for layername in layernamelist:
            ShapeFile.GetSpatialReferenceSystemType(worldsentinel2zoneshapefilepath, layername)
            geometrytype = ShapeFile.GetShapeFileGeometrytype(worldsentinel2zoneshapefilepath, layername)
            fielddictlist = ShapeFile.GetShapefileFieldList(worldsentinel2zoneshapefilepath, layername)
            ShapeFile.CreateEmptyShapeFile(newshapefilepath, layername, spatialreferencesystemwkt, geometrytype)
            ShapeFile.FieldDictListUpdateShapeFile(fielddictlist, newshapefilepath, layername)
            zonefeaturedictlist = ShapeFile.GetShapefileFeatureDictList(worldsentinel2zoneshapefilepath, layername)
            zonespatialreferencesystemwkt = ShapeFile.GetShapeFileSpatialreferencesystemWkt(
                worldsentinel2zoneshapefilepath, layername)
            newfeaturedictlist = ShapeFile.FeatureDictListTransformSpatialreferenceSystem(
                zonefeaturedictlist, zonespatialreferencesystemwkt, spatialreferencesystemwkt)
            ShapeFile.FeatureDictListUpdateShapeFile(newfeaturedictlist, newshapefilepath, layername)
        return

    @staticmethod
    def MultipleLayerCombination(
            sentinel2shapefilepath, countyshapefilepath, countysentinel2shalpefilepath, newlayername):
        """
        赛选某省县区的sentinel2图幅编号，生成相应的shape file
        :param sentinel2shapefilepath:
        :param countyshapefilepath:
        :param countysentinel2shalpefilepath:
        :param newlayername:
        :return:
        """
        sentinel2layercount = ShapeFile.GetShapeFileLayerCount(sentinel2shapefilepath)
        geometrytype = ShapeFile.GetShapeFileGeometrytype(countyshapefilepath)
        spatialreferencesystemwkt = ShapeFile.GetShapeFileSpatialreferencesystemWkt(countyshapefilepath)
        ShapeFile.CreateEmptyShapeFile(countysentinel2shalpefilepath, newlayername, spatialreferencesystemwkt, geometrytype)
        fielddictlist = 0
        countysentinel2featuredictlist = []
        for sentinel2layerindex in tqdm(range(sentinel2layercount)):
            sentinel2layername = ShapeFile.LayerIndexToLayerName(sentinel2shapefilepath, sentinel2layerindex)
            countylayername = ShapeFile.LayerIndexToLayerName(countyshapefilepath)
            intersectfeaturedictlist = ShapeFile.LayerIntersect(
                sentinel2shapefilepath, sentinel2layername, countyshapefilepath, countylayername)
            countysentinel2featuredictlist = countysentinel2featuredictlist + intersectfeaturedictlist
            if sentinel2layerindex == 0:
                fielddictlist = ShapeFile.GetShapefileFieldList(sentinel2shapefilepath, sentinel2layername)
        ShapeFile.FieldDictListUpdateShapeFile(fielddictlist, countysentinel2shalpefilepath, newlayername)
        ShapeFile.FeatureDictListUpdateShapeFile(countysentinel2shalpefilepath, countysentinel2shalpefilepath, newlayername)
        return

    @classmethod
    def GetSingleLayerFeatureDictListAccelerate(cls, args):
        print("start!")
        [shapefilepath, layername] = args
        singlezonefeaturedictlist = ShapeFile.GetShapefileFeatureDictList(shapefilepath, layername)
        cls.worldfeaturedictlist = cls.worldfeaturedictlist + singlezonefeaturedictlist
        print(len(cls.worldfeaturedictlist))
        print("end!")
        return

    @staticmethod
    def ShapeFileTransformSpatialreferencesystem(
            newspatialreferencesystemwkt, shapefilepath, layername, newshapefilepath, newlayername):
        """
        shape file 空间参考系转换（地理转投影）
        :param newspatialreferencesystemwkt: 
        :param shapefilepath: 
        :param layername: 
        :param newshapefilepath: 
        :param newlayername: 
        :return: 
        """
        shapespatialreferencesystemwkt = ShapeFile.GetShapeFileSpatialreferencesystemWkt(
            shapefilepath, layername)
        geometrytype = ShapeFile.GetShapeFileGeometrytype(shapefilepath, layername)
        ShapeFile.CreateEmptyShapeFile(
            newshapefilepath, newlayername, newspatialreferencesystemwkt, geometrytype)
        fielddictlist = ShapeFile.GetShapefileFieldList(shapefilepath, layername)
        featuredictlist = ShapeFile.GetShapefileFeatureDictList(shapefilepath, layername)
        ShapeFile.FieldDictListUpdateShapeFile(fielddictlist, newshapefilepath, newlayername)
        newfeaturedictlist = ShapeFile.FeatureDictListTransformSpatialreferenceSystem2(
            featuredictlist, shapespatialreferencesystemwkt, newspatialreferencesystemwkt)
        ShapeFile.FeatureDictListUpdateShapeFile(newfeaturedictlist, newshapefilepath, newlayername)
        return

    @staticmethod
    def WorldToChina(shapefilepath):
        """
        :param shapefilepath:
        :return:
        """
        layernamelist = ShapeFile.GetShapeFileLayerList(shapefilepath)
        zonenamelist = []
        for layername in layernamelist:
            layernamestrlist = layername.split("_")
            zonename = layernamestrlist[-1]
            zonenamelist.append(zonename)
            print(zonename)
        return

